Knock-on processes in superfluid vortex avalanches and 
pulsar glitch statistics 
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ABSTRACT 

A framework is presented for a statistical theory of neutron star glitches, motivated by 
the results emerging from recent Gross-Pitaevskii simulations of pinned, decelerating 
quantum condensates. It is shown that the observed glitch size distributions cannot 
be reproduced if superfluid vortices unpin independently via a Poisson process; the 
central limit theorem yields a narrow Gaussian for the size distribution, instead of the 
broad, power-law tail observed. This conclusion is not altered fundamentally when 
a range of pinning potentials is included, which leads to excavation of the potential 
distribution of occupied sites, vortex accumulation at strong pinning sites, and hence 
the occasional, abnormally large glitch. Knock-on processes are therefore needed to 
make the unpinning rate of a vortex conditional on the pinning state of its near 
and/or remote neighbours, so that the Gaussian size distributions resulting generically 
from the central limit theorem are avoided. At least two knock-on processes, nearest- 
neighbour proximity knock-on and remote acoustic knock-on, are clearly evident in the 
Gross-Pitaevskii simulation output. It is shown that scale-invariant (i.e. power-law) 
vortex avalanches occur when knock-on is included, provided that two specific relations 
hold between the temperature and spin-down torque. This fine tuning is unlikely in an 
astronomical setting, leaving the overall problem partly unsolved. A state-dependent 
Poisson formalism is presented which will form the basis of future studies in this area. 
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1 INTRODUCTION 

Pulsar glitches — the discrete, randomly-timed jum ps in the spin frequency of a pulsar 



law size and exponential waiting-time distributions (|Wong et al. 



200 ll : 



Melatos et al. 200; 



are characterised by power- 
These s tatistics point to a n 



underlying collective p rocess, similar to other complex systems such as neural and soci al networks ( Worrell et al.l |2002| ) , 
soil moisture b alance ( Porporato et al. 2004; Daly &: Porporato 200?!). electricity grids JCarreras et al 2002^). forest fires 



Turcotte*'l 999|). earthquak es ( Drossell 199d) and many mor e JGoss et al.lll989l : ISchonfisch fc de Rooj|l999l :l Cornforth et al. 



E 

l2005i : iBoerliist &: Hogewea 1991; Suki et aU 1994 : Wu et al. 201C). All of these systems respond to a global driver via local, 
stress-releasing interactions between discrete elements and may be modelled using many-body techniques. 

Pulsar glitches are caused by the unpi nning and near-simultan eous outward motion of many (up to ~ 10^*) quantised 
superfluid vortices in the pulsar inner crust ({Anderson fc Itohlll975l ). Individual, non-interacting vortices unpin stochasticall 



according to a Poiss on process, whose rate depends on the system temperature, local pinning forces IjHanggi et al.lll99l 



d; 



Hakonen et al. 1998h . local vortex density, and the relative velocity of the superfluid and pulsar crust. In addition to these 



factors, collective unpinning events are catalysed by unpinned vortices, which raise the unpinning rate of other vortices 
(nearby or distant, depending on the knock-on mechanism). These collective, avalanche-like events, in which vortices execute 
j erky, stick-slip motio n, proceed in a similar way to cascading failures in electrical grids and other large, networked systems 
( Carreras et al.|[2Q09l ). via an underlying threshold- driven branching process. They are characterised by a critical state in which 
the event size distribution is a power law with exponent —3/2. This simple power-law characterisation is most appropriate 
to a system in which the ev ent rate is constant. For systems like pulsars in which the event rate depends on some constantly 
evolving system parameter ( Dalv fc Porporaitol 2007 ). a more general formulation is required. 

The occurrence of a glitch brings the superfluid and crust closer to corotation and hence reduces the global unpinning rate. 
This feedback regulates the state-dependent vortex unpinning rate, which, when combined with the persistent electromagnetic 
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spin down of the su perfluid crust , leads to a self-organising system, whose evolution resembles that of a slowly driven sandpile 
or a forest fire (e.g. lOrossel 19961 ). That is, vortex unpinning decelerates the superfluid, reducing the crust-superfluid lag and 
hence temporarily decreasing the unpinning rate, which in turn increases over the long term in response to the electromagnetic 
spin-down torque. 

This paper asks whether a system of pinned vortices that unpin according to a variable-rate Poisson process, whose 
evolution is govern ed by pulsar-like rules, exhibit s unpinning avalanches like those observed. We build on the coherent noise 
model proposed by Melatos fc Warszawski ( 20091 ) . which successfully generates power- law-distributed glitches, with the aim 
of deducing the waiting-time distribu tion endogenously. instead of putting it in by hand. Vortex knock-on, as observed in 
quantum mechanical simulations (see Warszawski et al. 1 20121 )). provides a natural catalyst for unpinning avalanches. We 
propose a statistical model that self-consistently describes feedback between the superfluid and pulsar crust and introduces 
three novel features that reflect the latest understanding of pulsar interiors: stochastic unpinning, knock-on processes, and 
multiple pinning energies. The consequences of each feature are is assessed independently. We compare our analytic results 
to the output of Monte-Carlo simulations to test the underlying assumptions and explore regimes that cannot be treated 
analytically. 

The paper is organised as follows. In Sec. [2] we summarise the (un)pinning model of pulsar glitches and the challenges 
and paradoxes that it faces. Section [3] formulates vortex unpinning as a Poisson process. The feedback mechanism that relates 
vortex unpinning to changes in crust rotation is described in Sec. ID In Sec. [5] we summarise results from an asynchronous 
Monte-Carlo automaton, where the variable time step is governed by the shortest waiting time between vortex unpinning 
events, that realises the processes described in Sec. [2H11 We investigate the effect of pinning at one unique energy versus 
multiple energies in Sec. [S] and [7] respectively. In Sec.[H]we outline possible ways that unpinned vortices can catalyse unpinning 
avalanches, giving rise to a branching process, and describe how branching can be incorporated into both the analytic and 
Monte-Carlo treatments. It has proved useful to develop the Monte-Carlo and analytic descriptions side by side. In Sec. 19. II we 
describe how a master equation can be used to model the probability density function of the state of the system. The problem 
is for mulated in terms of a state-dependent process with Poisson transition statistics ( Dalv fc Porporatq 2007 : Wheatland 
20081 ). Conclusions and prospects for future work are contained in Sec. 1101 



2 VORTEX UNPINNING PARADIGM 



2.1 Standard model 

A pulsar's inner crust contains a neutron superfluid flowing through a crystalline lattice of nuclei. The superfluid rotates via 
the for mation of many (~ 10^^) quantised vortices, whose unperturbed configuration is an equally spaced, hexagonal Abrikosov 



lattice (|Donnellvlll99lh . Each vortex carries a quantum of circulation and generates a solenoidal velocity field. Macroscopically. 
the velocity fields from many vortices add vectorially to mimic rigi d body rotat ion. It is energetically favourable for a vortex 
core to sit atop a column of impurities (such as nuclear lattice sites ( JoneJl991b )). such that the volum e from which superfluid 
is excluded is minimised. The strength of th is phenomenon, known as pinning, is density dependent (jPizzochero et al.lll997l : 
Donati fc PizzocherollioO^ : |Pizzocherdl201ll ). I n addition, the extent of vortex pinn ing depends on the relative strength of the 
vortex- nucleus interaction and vortex tension I Link et al. 19931 : Haskell et al. 20121 ). 

As the pulsar spins down electromagnetically, pinning locks in the vortex lattice configuration, preventing the superfiuid 
from decelerating commensurately with the crust. Differential rotation between the ambient superfiuid and the pinned vortices, 
which necessarily corotate with the crust, exerts a transverse Magnus force on the vortex that, when sufficiently strong, unpins 
the vo rtex. The unpinning event and subsequen t outward vortex motion release angular momentum to the crust, causing a 
glitch (|Anderson fc Itoh|[l97i : lAlpar et al.lll98lh . 

The simplest form of the model outlined above is deterministic. That is, a vortex unpins when the shear between the 
container and superfiuid exceeds a critical threshold and stays pinned when the shear is subcritical. The Magnus force on 
pinned vortices grows linearly over time, and all vortices pinned at the same cylindrical radius experience the same force. 
Periodic glitches of equal size a re direct corollaries of this picture. Howe ver, periodic glitches of equal size are not seen in 
the data, except in two objects ( Melatos et al. 20081 : Estainoza et al. 2011 ) . The standard picture can be modified to include 
strong and weak pinning zones (jCheng et al.lll988l ) without treating local fiuctuations in pinning conditions at each vortex. 
However, whether or not a two-zone mean-field model is sufficient to reproduce the observed glitch size and waiting time 
distributions is an unanswered question. 



2.2 New physical ideas 

New physics needs to be added to the vortex unpinning paradigm in order to understand: (1) why some pulsars glitch and 
others do not; (2) why glitch sizes span up to four decades in an individual pulsar; (3) why there are large glitches involving the 
simultaneous unpinning of ~ 10^^ vortices, rather than a steady trickle of unpinning events; (4) why vortices skip over ~ 10^ 
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pinning sites before repinning; (5) if (how) large-scale inhomogeneities develop in a vortex lattice pinned to a homogeneous 
pinning array; and (6) how waiting times between glitches can be deduced endogenously from a model based on individual 
vortex unpinning. 

Our model addresses many of these questions (repinning is not specifically modellecQ) by incorporating the following 
three novel features: (1) individual vortices unpin explicitly according to a Poisson process whose rate is governed by the 
global superfluid-crust shear; (2) there is a range of available pinning energies, representing the myriad imperfections in the 
pulsar crustal lattice (we find that deep pinning wells are systematically overpopulated compared to shallow wells); and (3) 
unpinned vortices catalyse knock-on unpinnings at neighbouring pinning sites, which we model as a branching process. 

Domino-like knock-on unpinning, leading to collective vortex unpinning, is capable of explaining the observed scale- 
invariant dynamics of pulsar glitches. Num erical solutions of the G ross-Pitaevskii equation show that knock-on events arise 
natu rally in a pinned, rotating co ndensate ( Warszawski et al. 20121 '). and the resulting glitches span several decades in glitch 
size (| Warszawski fc Melatoj 120111 ). Knock-on events are triggered by vortex- vortex proximity, or acoustic radiation emitted 
during a nearby or rem ote repinning event. Knock-on stands in contrast to individual unpinnings characteristic of vortex 
creep IjAlpar et al.ll 19841 ) . although both types of events can occur simultaneously in the same system. 

The dynamics of individual vortices traversing an inhomogeneous pinning landscape is an active field of research. Gross- 



Pitaev skii s imulations show that v o rtices do not necessarily repin at the next available pinning site [see IWarszawski et al 
I 2012h and I Warszawski fc Melatod (120111)]. The exa ct trajectory depends on the superfluid-crust lag and the strength of 
pmmng. IjonesI (|l991al ). |jonesl l|l99a ) and lLiii^ (|2009l ) showed that there is a critical velocity below which a vortex travelling 
past a pinning site is immobilised. Its existence implies that the crust-superfluid lag necessary to overcome the pinning force 
must be large enough to propel a vortex past many pinning sites. Furthermore, in order to catalyse glitches of the size observed, 
vortices that unpin must move a distance comparable to the mean inter-vortex spacing (~ 1 mm) before repinning, which 

again involves bypassing many available pinning siteS; 

P revious models, such as the avalanche model in Warszawski fc Melatos 1 20081 ) and the coherent noise model in Melatos fc Warszawski 
I 2009[ ). suffer from two major drawbacks. The avalanche model requires inhomogeneities in the pinned vortex distribution on 
scales many orders of magn i tude greater than the average intervortex separation, contravening nuclear structure calculations 
i Donati fc Pizzocherol bood : lAvogadro et al.l 120071 ). The coherent noise model which is spatially uniform by construction, 
offers an ele gant way around this problem , by allowing pinning strengths to vary from site to site. However, in the form 
presented in lMelatos fc Warszawski l|2009l ). the waiting-time distribution is entered by hand, by fitting to observational data; 
it is not derived self-consistently. 



3 INDIVIDUAL VORTEX UNPINNING AS A POISSON PROCESS 

Pinned vortices move relative to the bulk superfiuid, experiencing a Magnus force per unit length equal to 

Fm = -pK X (vl - Vs) , (1) 

where k = {h/m)Cls is the quantum of circulation directed along the rotation axis f2s, p is the superfluid density, and vl — Vs 
is the velocity difference between the vortex line (which corotates with the crust) and the ambient superfluid. In the many- 
vortex limit, the superffuid mimics rigid-body rotation with angular velocity Q.^. Assuming that the vortices are straight and 
parallel and intersect the equatorial plane of the pulsar 0, we can write the Magnus force at radius R as 

Fm = pK7?AJ7cr-^ , (2) 
with Af2 = f2s — J^c = |vl — Vs|/i?, where fic is the angular velocity of the crust, a nd Aflcr is the m aximum differential 



angular velocity that a vortex pinned with energy Ep can withstand before unpinning IjAlpar et al.lll984l ). Aficr is a function 



^ How vortices repin remains an unsolved problem. One important question is how many pinning sites a vortex "skips over" before it 
repins. Upon examinin g many movies from Gross-Pitaevskii simulations of a pinned, decelerating condensate, whose results are reported 
in detail in IWarszawski fc Melatoj | |201 , we tend to observe that vortices skip over as many pinning sites as necessary to maintain 
roughly one Feynmann distance between each other. However, there are considerable fluctuations about this mean behaviour. Sometimes, 
a group of vortices unpin and leaves behind a gap in the Abrikosov array which is metastable and takes many avalanche "cycles" to 
fill. Other times, an unp inned vortex travels further than expected after grazing several pinning sites like a golf ball rimming the cup 
llWarszawski et al]|2012l ). None of these special situations are captured in the present paper and should be included in a more complete 
model. 

^ While it is traditional in neutron star literature to consider rectilinear vortices, and we continue that tradition here, simulations s how 
that at characteristic pulsar rotation rates, the flow in the inner crust is turbulent, and hence a tangle of vortices is likely fcr subota ct alJ 
I2OOOI : iPeralta et al]|2CI05l , I2OO6I : iMelatos fc PeraltallioOTi : lAndersson et al.ll2007l : iGlampedakis et aljliooj : iMelatos fc Peraltall20ia) . 
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of i5p, the superfluid coherence length, ^, and the pinning site separation 6, according to 



10~ 



Ep \ /I00fm\ nOfmX /10^'^kgm 



1 MeV J V 



10-* m 
R 



rad s 



(3) 
(4) 



Following Link et al. ( 1993h . we use Eq. ((2| to write the reduced pinning force, F' , on a pinned vortex that rotates differentially 
with the superfluid, as 

An 



F' = Fp- Fm = pnRAflc 



1 - 



(5) 



It is useful to describe vortex u npinning by ana logy with Kramers escape problem of a classical or quantum mechanical 
particle confined in a potential well ( Gardiner 20021 ) . That is, the Boltzmann escape or tunnelling probability increases with 
increasing system temperature and any local and global energy biases, including the differential angular velocity entering 
Eq. ((5]). Elence, we employ the Arrhenius formula to write down the unpinning rate, 6, of an individual vortex in a system 
with reciprocal temperature j3 — {ksT)'^, pinned with energy Ep, as (|Hanggi et al.lll990l : Ichevaliej|l993l) 

eiEp,'y)^roe-^^--'s-' , (6) 

where 7 = 1 — ASl/Aflcr parametrises the reduction in pinning energy due to differential rotation, and (r) — [9{Ep,'y)]''^ is 
the mean waiting time between unpinning events for a single vortex. The unpinning rate for a system of A^'v identical vortices 
is N.eiEp,^). 

The waiting times between unpinning events, r, are distributed according to an exponential probability density function 
(PDF), 



p(t; Ep, 7) = e{Ep,j) exp [-^(Bp, 7)r] 



(7) 



for a constant-rate (fixed 7) process. For a variable-rate process, S is a function of time through 7(t), and hence Eq. ((7)l 
becomes (IWheatland 



p(t; Ep, 7, t) = e[Ep,-y{t + r)] exp <^ - 



dt' e[Ep,j{t + t')] 



(8) 



The integral in the exponent in Eq. ([8]) accounts for the increasing likelihood that unpinning occurs as time passes due to the 
accumulating crust-superfluid lag. 

In the simple case of uncorrelated unpinning events, the PDF of event sizes is trivially 

p{n) = S{n - 1) , (9) 

where n is the number of vortices that unpin. Equation (|9]) highlights the fact that there are no avalanches. Individual 
unpinning events can still cluster in time, as in any Poisson process. Unlike knock-on-mediated avalanches (see Sec. [8]), 
however, where even a system-spanning avalanche lasts much less than the typical time between events, and it is obvious 
what constitutes a single avalanche, there is no preferred time-scale for identifying complete events (except the trivial one of 
recording unpinnings one by one, as above). In a pulsar, many individual unpinnings occur in a typical observational time 
window At. Binning events according to the observational resolution At, we obtain the familiar Poisson result 

(dAt)" 



p{n) = exp(-6'At)- 



(10) 



which tends to a Gaussian distribution for large At, and 9 remai ns the unpi nning rate for a single vortex. The Gaussian 
conflicts with data in all pulsars except Vela and PSR J0537-6910 IjMelatos et al.,,2008, '). 



4 FEEDBACK AND SELF-REGULATION 

A glitch brings the coupled superfluid- crust system closer to corotation. Immediately after the event, the superfluid decelerates 
and the crust accelerates, reducing the vortex unpinning rate by increasing 7 [see Eq. (|8])]. After a longer delay, the electro- 
magnetic torque decelerates the crust (without affecting the superfluid), gradually reducing 7 and hence increasing 6. Spin- up 
of the crust in response to downward steps in superfluid rotation ensures that the glitch mechanism is self-renewing. That 
is, each time a glitch occurs, the lag between the crust and superfluid is partially reset (almost never to zerojf], reinitiating 
the accumulation of lag that leads to the glitch. Without feedback, the lag would continue to grow, eventually rendering 

^ Since no correlation is observed between glitch size and waiting time, it can be deduced that the lag between the superfluid and crust 
rotation is not reduced to zero by a glitch. 
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Figure 1. Results from a mock simulation, in which power-law distributed glitch sizes and exponential waiting times are superposed on 
a uniformly decelerating crust. Equation Ulll l describes how changes in f2c and fig are related. Left: Angular velocity as a function of time 
HcC*) for the crust (Hci black curve) and superfluid (Qg, grey curve; arbitrarily shifted down by 0.706 for comparison). The superfluid 
spins down in a step-wise fashion, corresponding to vortex unpinning, whilst the crust angular velocity decreases linearly with time, 
accelerating instantaneously when the superfluid decelerates. Right: Differential rotation parameter as a function of time, 7(t), for the 
same interval graphed in the left panel. ■y{t) decreases linearly between unpinning events, and increases in discrete jumps when vortices 
unpin. Simulation parameters: Is = Ic- 
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Figure 2. Results of a mock simulation showing how the PDF of observed glitch sizes depends on the observation time window for 
binning. At. In the underlying simulation, a power-law distribution of glitches (index —3/2), separated by exponentially distributed 
intervals, is superposed on a background of constant deceleration. Left: Angular velocity of crust as a function of time Qc{t)- Right: 
PDFs of 'observed' glitch sizes for four time windows. At = lO"'^'^^, 1D~^'^^, 10"'^ '"' and 10""'^^ {black to light grey curves respectively). 
The glitch-finding algorithm is described in Sec. 15.21 



pinning impossible. On very long time-scales (> 10^ yr), therefore, the system is not stationary. However, the slow electro- 
magnetic spin-down rate, dQc/dt (~ 10~^^ Hzs~^, for a typical pulsar) suggests that, over the forty years of historical pulsar 
observations we can treat Qc as constant. 

The crust responds instantaneously to changes in the superfluid angular velocity according to 

^^^ = -^=^+^- (11) 
where Ic and h are the m oment of inertia of the pulsar crust [and the charged proton and electron fluids that are tightly 
coupled to it l Alpar et al.l 1984)] and the superfluid (notionally, if it were a rigid body of the same density) respectively, and 
A'^c is the electromagnetic torque on the pulsar crust Q In practice, the time scale on which changes in the superfluid angular 



^ We remind the reader that the observed deceleration of a pulsar is a compromise between A'^c and the time-averaged spin-up torque 
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Figure 3. Monte-Carlo simulations of a vortex system with a unique pinning energy, torque feedback described by Eq. and no 

knock-on avalanches. The automaton uses the rules described in Sec. (5] for a system of A^v = 10'^ vortices with Ep = Eq = 1.0, 
N^/I^ = 10-° l, a = 10-'', /s//c = 1.0 and Nt = 10*^. Results are reported as time series and PDFs. The PDFs are constructed by 
sampling the time-varying quantity, and making a histogram of the sampled values. The glitch-finding algorithm, used to construct 
p(Af!), is described in Sec. 15.21 Top left: Angular velocity as a function of time, f2c(t), for the time interval 2.519 <t< 2.557 (solid 
curve). The dotted curve is a linear fit to Clc{t) with slope —0.390. Top right: PDF of jumps AQc for three glitch-finding time windows. 
At = 10-2 0^ 10-1-5 and lO-^". Centre left: Time series of waiting times since previous glitch, t, versus glitch epoch t. Centre right: 
PDF of T on log-linear axes {solid curve). The dotted curve is an exponential least-squares fit with exponent 9 = —0.42. Bottom left: 
Shear parameter as a function of time, 7(<). Bottom right: PDF of 7. The dashed vertical line indicates the position of 7oq = 0.923, 
which is offset from the mean {7) = 0.922. 
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momentu m are communicated to the crust is governed by the crossing time of a Kelvin wave along the length of a vortex 
[~ 10~^ s I Epstein fc Bavm 1992l l]. We can rewrite Eq. Ulip in terms of changes to the global shear parameter 7, 



1 



- 1 + 



+ 



Ic 



(12) 



It should be noted, that interactions between the viscous component of the internal fluid (composed predominantly of protons 
and electrons) and the neutron superfluid also result in a coupling to the crust, hence influencing the speed and strength of 
coupling. The likely time scale of suc h an interacti on is the Kelvin wave crossing time of the star, which is also equivalent 
to the mutual fri ction coupling time |Mendel] 1991), which has been studied recently in a hydrodynamical context in glitch 



recovery models (|van Evsden fc Melatos 



201' 



1 



Strictly speaking, a constant-rate {dy/dt = 0, constant /?, Ep and Fq) Poisson process does not capture the self-organising 
behaviour that is essential to regulating the crust-superfluid lag, since changes in the lag are equivalent to changes in the 
unpinning rate. Fluctuations in 7, caused by glitches, render vortex unpinning an inherently non-constant-rate process. The 
absence of a reservoir effect in glitch data (glitch sizes are not correlated with waiting times) indicates that glitches never 



comp letely nullify the lag, so the stress released in each glitch is a combination of recent and historical build up (|Melatos et al 
2003). In Sec. 19.11 we discuss an analytic formalism for describing a variable rate, or state- dependent, Poisson process. For 



now, however, we investigate the behaviour of a system with a constant rate, employing Eq. (Illll to ensure that unpinning 
balances the spin down of the crust on average over the long term. 

The change in angular momentum of the superfluid when one vortex unpins and moves outward a fraction a of the stellar 
radius is given by (see Appendix IX)) 



(13) 



Gross-Pitaevskii simulations suggest that ctR approximately equals the mean inter- vortex spacing IjWarszawski et al.ll2012l ) 
with 



1 / K \ i/2 

'r vm) 



(14) 



whether the number of vortices exceeds or is exceeded by the number of pinning sites IjWarszawski fc Melatosll201ll '). In order 
to find the asymptotic spin-down rate, dQs/dt, we now consider a constant-rate process, with 9[Ep,'y{t)] defined in Eq. Q. 
That is, for d"f/dt = we get 



iVc//c 



dO,s 

~dr 1 + Is/Ic 



(15) 



Although in practice dy/dt is non-zero, self-regulation, via Eq. (|lip ensures that the unpinning rate results in the commensurate 
spin down of the superfluid with the crust. On average, the superfluid spins down at a rate proportional to the mean vortex 
unpinning rate, viz. 

^ = ^iV.roexp(-/3i5p7oq) • 

Combining Eq. (|15|) and (|16|) . and solving for 7, we have an estimate for the steady-state shear parameter, 7eq, 



(16) 



7cq 



■In 



/s//c 



PEp VToA^vAL 1 + J,/Jc 



(17) 



which ensures that superfluid spin down by unpinning matches the rate dictated by the electromagnetic torque. 

The left panel of Fig. [T] illustrates how ilc and Sis change when allowed to evolve according to Eq. pip for a set of 
artificially generated, power-law-distributed glitch sizes, imposed at exponentially distributed intervals on a background of 
constant deceleration. Between unpinning events. Sis (grey curve) does not change, whilst Qc {black curve) decreases linearly 
with time at a rate Nc/Ic- When vortices unpin, the Sis curve steps down instantaneously, matched by an instantaneous 
acceleration of the crust. In the right panel we graph 7(t) over the same time interval; between unpinning events, 7(4) 
decreases linearly at a rate Afc/(/cASlcr). It increases instantaneously by (ASl — ASls)/ASlcr when a vortex unpins, with 
ASl = AL//c. 



exerted by the decelerating superfluid, which is composed of both smooth and discrete jumps. Therefore, we cannot merely write 
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Figure 4. Monte-Carlo simulations for the vortex system in Fig.[3l but with a range of available pinning energies Ep S [Eq — Ai?p, Eq + 
Ai?p]. The automaton rules are described in Sec. [5] Results are reported as time series and PDFs. Top left: Angular velocity as a function 
of time, Qc{t), for the time interval 1.084 < t < 1.158 (solid curve). The dotted curve is a linear fit to Qc{t) with slope —0.390. Top 
right: PDF of jumps, AS7c- Centre left: Time series of time intervals between unpinning events, t. Centre right: PDF of r on log-linear 
axes (solid curve). The dotted curve is an exponential least-squares fit with rate parameter 8 = —0.42. Bottom left: Shear parameter as 
a function of time "((t). Bottom right: PDF of 7. 7cq = 0.751 falls outside the plotted range; it is larger than the average (7) = 0.697. 
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5 MONTE-CARLO SIMULATIONS 

In the remainder of the paper, we implement the physics in Sec. [3] and |4] in an automaton to study the open questions described 
in § II. We vary the rules governing the vortex unpinning rate and the distribution of pinning energies and analyse the output 
statistically, aided by a glitch-finding algorithm, which compiles glitch waiting-time and size distributions (see Sec. 15. 2[) . 

The automaton output and the analytic estimates given in Sec. [4] are compared against 7cq, given in Eq. (|17[l . and the 
mean of the simulated 7 distribution. Equation (|17|l only applies for uniform pinning energies, Ep = Eq. Later, in Sec. 19.11 
we write down a master equation [Eq. (|37|l] that corresponds to the uniform-pinning scenario. We extend the simulations and 
theory to encompass a range of Ep in the next section. 



5.1 Asynchronous automaton 

In an asynchronous cellular automaton, the time step is a stochastic variable that depends on the state of the system 0. 
Each step in the automaton is governed by the following rules. 

(1) For each vortex i, a time until the next unpinning is drawn from p{T;Ep,t) [defined in Eq. ((8|], with the initial 
condition "f{t — 0) =70. 



1999 ) 



2) The system is advanced in time by Tmin = min{ri}, and the global time counter becomes t„ — tn-i + Tmin (jBiane et al 



(3) Qs is decremented by Ails = AL/Is, and fie changes by AQc ~ —AL/ Ic + TmhiNc/ Ic [see Eq. 1)11^ ]. The shear parameter 
becomes y{t + Tmin) = jit) + (Afic — AQs)/ AQcr- 

(4) The unpinned vortex is assigned a new waiting time (and a new pinning strength if AEp 7^ 0), drawn from p(r; Ep, t). 

(5) For all vortices that did not unpin in the previous interval, the time until the next unpinning event becomes 

{n - Tmin) exp [-l3Eo-i{t + Tmin)] / exp [-PEo-fit)] . (18) 

(6) Repeat steps 2-5. 

Rules 4 and 5 define what is known in the literature as a time-driven asynchronous automaton; instead of updating every 
cell at each time step, cells are updated in the order that Poisson events affecting them actually occur. The stationary states 
of an automaton do not depend on whether the update rules are synchronou s or asynchronous, but this is not true of their 
basins of attraction (iHuberman fc Glancel 19931 : Schonfisch fc de Roos 1999h . Asynchronous update rules are both efficient 



and fiexible. When the event rate is slow, an asynchronous algorithm avoids unnecessarily small time steps, while the lack of 
a minimum time step means there is no upper bound on event rates. 

Rule 5 is necessary because each event changes the state of the system and hence the unpinning rate for all vortices. 
Equation (|18p encodes the response of the unpinning rate to the new state of the system, by stretching the remaining fraction 
of the waiting time according to the ratio of new and old rates. That is, the vortex with the shortest waiting time unpins after 
Tmin time units. Another vortex j would have unpinned in r, — Tmin later without any adjustment. However for the remainder 
of the waiting time, the unpinning rate changes and, to compensate, we stretch Tj — Tmin by the ratio of the new and old 
unpinning rates, Fo exp[— /?i5o7(t -I- Tmin)] and Fo exp[— /?i5o7(t)] respectively. 

An alternative approach would be to reset the waiting times of all vortices at the beginning of each new time interval. 
Such an approach agrees with the spirit of a Poisson process, in which time intervals of any length are independent, but it is 
not favoured by other authors. We follow the majority custom here. 

In all forthcoming plots, quantities are given in arbitrary units. 



5.2 Glitch classification 

A robust, automated glitch-finding algorithm is essential for compiling and studying glitch statistics. The simplest approach 
is to calculate the accumulated spin up, Afic, during some time window of width At, and then construct the probability 
density function (PDF) of Afic, p(Af2c), by sliding this window along the time series Q.c{t). 

The shape of p[AQ.c) depends strongly on At. Figure [2] graphs f2c(i) for a mock system, in which glitches with power- 
law-distributed sizes (exponent a — —3/2) and exponentially distributed waiting times are superposed on linear deceleration, 
viz., 

n4t) = n^o) + n,t + J2 ^^sii) ' (i9) 

tetg 



^ The 'best choice' for the order in which cells are update d depends on the spec ific application; it has been shown that the choice of 
update scheme can greatly influence the simulation output l lCornforth et aUl2005l ). 
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where AS7g(t) represents the intermittent glitches and is non-zero only when a glitch occurs at f = tg (in arbitrary units). 
The left panel graphs n^t). The right panel graphs p(Af2c) for At = 10~^■■^^ lO"^'^^, 10~" '^\ 10"°-^^ {dark to light curves 
respectively). For the two smallest values of At, p{AQc) is well represented by a power law, but the best-fit power-law index 
(a = —1.31 for At = 10~^'^*) does not match exactly the distribution from which the glitches were drawn (a = —3/2). For 
larger values of At, p(Ar2c) peaks about a mean value. This behaviour is a direct result of the central limit theorem. When 
the time window is much larger than the mean inter-glitch waiting time, each window brackets several glitches, and hence 
Afic, which is the sum of several draws from the same underlying size PDF, populates a Gaussian distribution. Therefore, we 
conclude that a narrow window with At < 9^^ best characterises the underlying glitch behaviour. 



6 PINNING AT A UNIQUE ENERGY 

The simplest possible form of the vortex (un)pinning theory of glitches treats the pinning strength as uniform throughout the 
pulsar. Naively, this scenario should produce quasi-periodic glitches with a narrow range of sizes, since each vortex sustains 
roughly the same level of stress before unpinning. In this section we test this idea using the automaton described in Sec. [5] 

Results from Monte-Carlo simulations for single-energy pinning {Eo — 1, AEp — 0) are shown in Fig.O The electromag- 
netic torque is A'^c = 10~'''^/c, each unpinned vortex travels a distance aR = 10~*ii, the crust and superfiuid moments of 
inertia are taken to be equal, (7s//c = 1), and there are = iC vortices. The rules are enacted A't = 10^ times. The graphs 
show the angular velocity of the crust as a function of time ilc{t), and its PDF, ^(Ailc), waiting times between unpinning 
events r as a function of time, and their PDF p(r), 7 as a function time, and its PDF ^(7) (from top to bottom, left to right). 
The PDFs are constructed by sampling the quantity at each time step and constructing a histogram of the sampled values. 
Calculation of ^(ASlc) is discussed in Sec. 15.21 

The time series of Sic exhibits restless behaviour, including spin-up events. The dotted curve is a straight line with slope 
—Nc/{2Ic), which is the expected spin-down rate for 7 = 7eq and /a/-^c ~ 1. Although the simulated Sic wanders away from 
the grey curve, —Nc/{2Ic) remains a good approximation for the global spin-down rate over long time periods. The size 
distribution, ^(ASlc), is strongly peaked, with mean and standard deviation 0.40 and 0.02 respectively. It does not span many 
decades, unlike real observed glitches. Further evidence that fluctuations in Q,c due to glitches do not push it far from the 
time-averaged trend is found in the ^(7) curve, for which 7cq = 0.923 [dashed vertical line, taken from Eq. (|40|l ] lies within 
the narrow range of values (mean and standard deviation 0.92 and 0.0023 respectively) visited by the automaton. However, 
we note that 7cq does not coincide exactly with the mode of ^(7). 

The narrow distribution of 7 values leads naturally to an almost exponential distribution of waiting times because the 
activity approaches a constant-rate Poisson process, with mean rate given by S ~ A'^vFq exp(— /37eqi?o) {grey curve in middle 
left panels of Fig. [3] and |4]). The dotted curve is an exponential least squares fit, whose slope differs from — /37cq£o by 0.053. 

We conclude that a state-independent Poisson process captures well the behaviour of the single-pinning-energy system. 
However, such a system does not give rise to a broad distribution of event sizes, like those observed astronomically. In the 
following sections, we repeat this analysis for systems with (i) a wide range of pinning energies, and (ii) collective triggers for 
unpinning avalanches (Sec. [7] and |8] respectively). 



7 PINNING AT MULTIPLE ENERGIES 



The strength of pinning within a pulsar is likely to be heterogeneous. Small-scale (microscopic) het erogeneities arise from 
defects in the crystalline lattice of the crust, in the f orm of either vaca ncies, dislocations, or impurities i Pizzochero et aLlll997 : 
Blasio fc Lazzari 19981 : Donati fc Pizzochero 20031 : Pizzochero 2008h . Large-scale (macroscopic) fluctuati ons in the pinning 
strength may resu lt from grain boundaries in the crust, creating reservoirs of strongly pinned vortices ijAlpar et al.lll984l : 
Cheng et al. 1988h . We approximate the effect of small-scale variations by choosing the pinning energy at each pinning site 
from a distribution function (l){Ep), defined as 



0(£p) = {2Eo)-'H{Eo - AEp - Ep)H{Eo + AEp - Ep) 



(20) 



where -Eo is the mean pinning energy, AEp is the half-width of the distribution, and H{-) is the Heaviside step function. 
Whenever a vortex repins, its new pinning energy is chosen at random from (l){Ep). 

Heterogeneous pinning intro duces a new statistical quantity, namely the distr ibution of occupied pinning energies, g{Ep, t). 
This is not the same as 4>{Ep) ( Newman fc Sneppen 19961 : Melatos et al. 20081 ): shallow pinning wells are less likely to be 
occupied than deep ones, in a statistical sense. In Sec. 17.21 we estimate g{Ep,t) for a stationary system analytically. 
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Figure 5. Schematic of tiow vortices accumulate over time at strong pinning sites, biasing the distribution of occupied pinning sites, 
g(Ep), towards large Ep. Since weakly pinned vortices repin more frequently than strongly pinned vortices, and a vortex repins with a 
strength chosen at random from a uniform distribution, (f){Ep), over time, weak pinning sites are depopulated in favour of strong pinning 
sites. The panels depict this process taking place over four episodes, starting with uniformly-distributed, random pinning strengths {top 
left), and ending with a preference for occupying strong pinning sites {bottom right). In each panel, the vortices which are about to unpin 
(and repin uniformly in the next panel) are coloured grey. 




Figure 6. Pinning energy distribution of occupied sites, g{Ep), for the Monte-Carlo simulation reported in Fig.|4](soK(i curve), overplotted 
with gcq{Ep) from Eq. I I26II {dashed curve). The solid grey curve is g{Ep), calculated using 7 = (7), taken from p{-f) {bottom right panel 
of Fig.|4]l. The light grey histograms display the instantaneous g{Ep,t), at five equally spaced times. 



7.1 Automaton output 

Figure |4] graphs results from Monte Carlo simulations for heterogeneous pinning {AEp/Eo = 0.5, Eq — 1, Is/Ic = 1, 
Nc/Ic = lO""'^, a = 10"*, A'v = 10*, Nt = 10^). The graphs show the angular velocity of the crust as a function of time 
r2c(t), the PDF of spin-up events, p(Ar2c), waiting times between unpinning events r as a function of glitch epoch, and the 
PDF p(r), and 7 as a function of time and the PDF (from top to bottom, left to right). 

Qualitatively, the spin-down curve {top left panel) does not differ greatly from the unique pinning energy case {top left 
panel of Fig. [3]). The PDF of Af^c confirms that heterogeneous pinning strengths are insufficient to significantly broaden {i.e. 
across several decades) the Af^c distribution; the range of AQ,c is less than one decade (10"" '' < Afic < 10°'^), slightly but 
not much wider than ^(ASlc ) in Fig. El 

The validity of a state-independent (7 — 7cq) description of the system is confirmed by the relatively narrow distribution 
of 7 (0.79 < 7 < 0.8). The waiting times (middle right panel) are well represented by an exponential distribution with mean 
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rate given by the equilibrium unpinning rate, 6'cq = 0{Ep, jcq) (solid grey curve), which in turn closely matches the exponential 
fit plotted as a dotted curve. 



7.2 Excavation 



Vortices unpin readily from weak pinning sites, and with difficulty from strong pinning sites, as illustrated schematically in 
Fig.[5](the sequence is from A — D). This tendency 'excavates' the PDF of occupied energies so that the majority of vortices sit 
on strong pinning sites. In the initial state, immediately following a system-spanning glitch, vortices are pinned at sites with 
uniformly distributed energies (as in frame A). Vortices that subsequently unpin are coloured grey. Weakly pinned vortices 
are more likely to unpin and repin at randomly distributed energies (as in frame B). A preference for strongly pinned vortices 
emerges (frame D). This behaviour, to which broadly di stributed glitch sizes were attributed, was also observed in the coherent 
noise model described in Melatos fc Warszawskil (|2009l ). 

Let g{Ep, t)dEp be the fraction of vortices at time t pinned at sites with pinning energies in the range {Ep, Ep + dEp). In 
a time dt, in the energy range (Ep, Ep + dEp), a number of vortices equal to 6(Ep,"f)g(Ep, t)dt unpin. Some repin in the same 
range (Ep,Ep + dEp), in numbers proportional to the fraction (j>(Ep)dE of available sites in that range, and some repin at 



other sites with energies outside this range. Inspired by models of atomic hopping i n glasses llSouchaud et al1ll995 



Rinn et al. 



2001), this stochastic process can be described statistically via a master equation I Monthus fc BouchaudI 19961 : HeadI 2000l ) 

dg(Ep,t) 



dt 



= -e[Ep,-/(t)]g(Ep, t) + uj(t)(t>(Ep) 



(21) 



where ij{t) is the unpinning rate at time t averaged over all occupied sites, 

uj(t)= dE'pe[E'p,^(t)]g(E'p,t) . (22) 
Jo 

Equation (|2ip is solved supplemented by the initial condition g(Ep, 0) = (j>(Ep) ( Bak fc Sneppen 19931 ) by taking the Laplace 
transform (see Appendix|B]). We kn ow from observation s — e.g. the universal exponential waiting-time distributions observed 
in pulsars on the time-scale of years ijMelatos et al — that the glitch statistics and therefore presumably the underlying 

g(Ep,t) reach stationarity fairly quickly. From Eq. (|21[1 . in the steady state [dg(Ep,t) / dt = 0], u)(t) must tend to a constant 

l^eq(/3,7cq), wlth 



iVvFo 



2A£p/37eqe 



dE' e"'"'^''' <t,(E') 



(23) 



(24) 



where Eq. H24p is obtained from Eq. (|23|) by substituting Eq. (|20[) . Hence the stationary distribution gcq(Ep) — g(Ep,oo 
takes the form 

-HEp) 



fleq(-Ep) = 



iVvFo 



o^7cq-Ep 



(25) 
(26) 



Note that the stationary state always exists provided that the integral in Eq. (|23|l is finite. 

In practice, g(Ep,t) is inherently time dependent, because 7 is time dependent. Figure [6] graphs g(Ep,t) from the Monte- 
Carlo simulations discussed in the previous section (AEp/Eo = 0.75) at five equally spaced times (solid, light grey curves). 
We construct the stationary gcq(Ep) (solid black curve), using 7 — (7) taken from p('y) (bottom right panel of Fig. |4)), by 
averaging g(Ep) from 50 snapshots in time. For comparison, we also plot 5cq(i5p) from Eq. (|26p (dashed curve). The light grey 
curves display the instantaneous g(Ep,t), at five equally spaced times. For both the average and instantaneous (?cq(£p), there 
is a systematic overpopulation of strong pinning sites compared to weak pinning sites, illustrating the excavation property 
illustrated in Fig. [5] We attribute the discrepancy between g(Ep) and Eq. (|26p to the fact that Eq. (|21|l does not allow for 7 
to vary in response to glitches. We also note that, in results not presented here, the simulated geq(Ep) is better approximated 
by Eq. (|26p when synchronous automaton rules are used. 



8 VORTEX AVALANCHES 

In Sec. [Hill we find that when each vortex unpins according to an independent Poisson process and inter- vortex correlations 
are ignored, vortex unpinning avalanches with a broad range of sizes do not occur, whether or not multiple pinning energies 
are available. Indeed, vortices unpin one by one and the only way to get events that involve n > 1 vortices is to sum them 
together over an observational window At (see Sec. l5.2|) . Therefore, to explain the catastrophic, near- instantaneous unpinning 
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Figure 7. Results from Monte-Carlo simulations of a system with asynchronous update rules and knock-on between vortices. Results 
are reported as time series and probability density functions (PDFs). Top left: Angular velocity as a function of time, Clc{t), for the time 
interval 130.7 < t < 134.6 {solid curve). The dotted curve is an empirical linear fit to Clc(t) with slope —0.40. Top right: PDF of number 

of vortices that unpin in each avalanche {solid black curve); a power-law fit with power-law index —1.98 is overplottcd {dotted curve). 
Bottom left: Time series of 7. Bottom right: PDF of 7. Simulation parameters: A^v = 500 vortices with pinning energies in the range 
Ep e [Eo - AEp,Eo + AEp], Nc/h = lO""-!, a = IQ-^, h/h = 1.0, Eq = 20.0, AEp/Eo = 1 and = 5 x 10^. 
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Figure 8. Left: PDF of number of vortices that unpin in each avalanche, n, for the same simulations as described in Fig.O This should 
be compared to the top right panel of Fig. [T] where glitch size is given as AQ. Right: g{Ep) for simulations reported in Fig. ^ {solid 
histogram). The bias towards weak pinning sites contradicts the excavation concept described in Sec. 17.21 



3.0 - 
2.8 - 
2.6 - 
o 2.4 - 
2.2 - 
2.0 - 
1 .8 - 



4.0 4.5 5.0 5.5 6.0 6.5 



Figure 9. Best-fit power-law indices for the glitch size (measured in number of vortices unpinned) distribution resulting from Monte 
Carlo simulations which include knock on between unpinned and still-pinned vortices. Knock on is parametrised by a boost in the 
unpinning rate due to unpinned vortices moving through the pinned vortex lattice. The PDF is constructed as a function of this boost, 
A7/A70 for Monte Carlo simulations; A70 is given in Eq. II28I I. Simulation parameters: A^v = 500 vortices with pinning energies in the 
range Ep S [Eq - AEp, Eq + AEp], N^/h = lO-o i, a = lO"*, I^/I^ = 1.0, Eq = 20.0, AEp/Eo = 1 and A^t = 5 x 10^. 



of vortices in groups that is responsible for pulsar glitches, collective unpinning mechanisms must be considered. In this 
section we quantify how unpinning events multiply when a recently unpinned vortex triggers unpinning of another pinned 
vortex, laying the basis for a domino-like knock-on process. Results from first-principles, quantum mechanical simulations of 
single vortex unpinning and repinning (Warszawski ct al. 2012) motivate this study. They point to two knock-on mechanisms: 
unpinning by the Magnus-like proximity efTect, and unpinning by acoustic waves. In what follows, we focus on the former. 

According to the Feynman-Onsager quantisation condition, the change in the superfluid speed at one of the vortices. 
An, when the sepa ration between a moving vortex and its nearest neighbour decreases from dp to edp (e < 1) is given by 



I Cheng et allll988l ) 



Just as in the global shear lowers the pinning barrier via 7, Eqs. (|27[) and ([1]) demonstrate that a local over-density of vortices 
(relative to the unstressed lattice configuration) lowers the pinning barrier too. Locally the shear parameter 7 changes from 
7 to 7 — A70, with 

K 1 — e 

= 27rdpi?AJ7„ . • ^28) 

Suppose, for the sake of illustration, that all the vortices are evenly spaced by dp and one vortex unpins; call this a 
trigger event. The unpinned vortex moves radially, approaches a nearby vortex, and reduces the latter's pinning barrier. The 
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time-scale for this to occur is Tk = dp/uv ~ 10 s, where = RAQ is the radial component of the unpinned vortex velocity. 
During this time interval the unpinning threshold is reduced by A7, such that the unpinning rate for a single vortex becomes 

A(£;p) = roe-''^p(^-'^^' . (29) 

So, each trigger unpins, on average, A(i5p)rk vortices. In practice, to see whether or not this formulation of knock on results 
in unpinning avalanches, Eq. (|29|1 should be incorporated into the automaton rules described in Sec. 15.11 which we do in the 
next section. 



8.1 Automaton rules 

In this section, we modify the automaton rules described in Sec. [5] to include knock-on. 

(1) For each vortex i, a time Ti until the next unpinning is drawn from p{Ti,Ep, t), with the initial condition "/{t = 0) = 70. 

(2) The system is advanced in time by the minimum waiting time Tmin ~ min{ri}, t„ — t„-i + Tmin. 

(3) Increment Qc by Afic = TminA^c/Zc. The state parameter 7(t -I- Tmin) becomes ■y{t + Tmin) = 'y{t) + (Afic — AS7s)/Af2cr. 
Now allow the avalanche to proceed. 

(a) Ask each vortex (other than the original trigger) if it unpins, with probability 1 — exp(— Ark) [the probability of one 
or more events during an interval rk, characteristic of the knock-on process, in a Poisson process]. A7 is parametrised in 
terms of A70 [Eq. ((28)) ]. 

(b) If at least one extra vortex unpins, ask all the still-pinned vortices if they unpin until none do and the avalanche is 
over. The avalanche size is proportional to the total number of vortices that have unpinned, n. 

(4) Decrement fls by Afls = nAL/Is, so that Sic changes by AQc ~ —nAL/Is- The shear parameter 7(t-|-rmin) is updated 

to 7(t + Tmin) = 7(i) + (Afic - AQ,s)/AVL^r. 

(5) All unpinned vortices are assigned new waiting times (and new pinning strengths if AE^ 7^ 0), drawn from piji, Ep,t-\- 
Tmin). Notc that this is the second 7 update within a single time step. 

(6) For all vortices that did not unpin in the previous interval, the time until the next unpinning event for the i^^ vortex 
at the n*'' iteration, Ti,„ becomes 

(t, - Tn^in) exp [-/3£o7(i + -Tmin)] / cxp [-/?£'o7(t)] , (30) 

which stretches the remaining fraction of the waiting time for still-pinned vortices in response to the new state of the system. 

(7) Repeat steps 2-5. 

The above rules differ from those given in Sec. 15.11 Each unpinning is treated as a trigger for subsequent unpinnings, which 
occur at a higher mean rate, causing an unpinning avalanche. In Sec. 15.11 by contrast, the subsequent unpinnings are not 
included, such that n = 1 always. 



8.2 Output statistics 

The impact of knock-on on the event size distribution depends on the relative strength of pinning and the knock-on variable, 
A7 [we discuss the physical implications of different A7 in Sec. (|9.ip ]. To understand why, first consider a weak, single- 
pinning-energy scenario, in which a relatively small AO. (high 7cq) is sufficient to unpin vortices. In this case, A7 is small 
compared to 7eq, and hence the increase in the unpinning rate due to the knock-on effect is small. On the other hand, for 
strong pinning, 7eq is smaller, rendering knock-on more effective. In this way, either increasing _Ep or decreasing A7 drives 
the system towards avalanches. 

Results from Monte-Carlo simulations that employ the asynchronous update rules described above (Sec. 18.1)) . are presented 
in Fig. [7| In the top left, we graph Q.c(t) [solid black curve), overplotted with a linear fit [dotted black curve), as well as the 
curve that corresponds to uniform spin down [slope —Nc/[2Ic), solid grey curve]. Unlike the corresponding plots in Fig. [3] 
and m Slc(i) is punctuated by abrupt jumps, caused by near-instantaneous knock-on events (the characteristic knock-on 
time scale, quantified in Sec. 18.31 satisfies Tk <C (t)), followed by periods of uniform deceleration. On long time scales, the 
spin down matches that expected for a container filled with an unpinned, smoothly decelerating superfluid [Eq. (|lip with 
d^lc/dt — d^ls/dt and Is/Ic ~ 1, gives dQc/dt — Nc/[2Ic)], as shown by the agreement between the fitted [dot-dashed) and 
linear [grey) spin-down curves. 

The top right panel of Fig. [Tjgraphs p(Af2c), calculated using three At fSec. 15.2]) . Interestingly, p[AQ,c) spans over three 
decades in AQ,c, as compared to the narrow AQ,c range without knock-on (see Fig. [3] and [4]). There is a plateau between 
0.1 < AQ,c < 1.0 for the two smaller values of At. Further investigation is required to understand this feature. 

We also measure glitch size by the number of vortices that unpin during an avalanche, n. Quantitative agreement between 
this approach and the glitch- finding algorithm in Sec. I5.2l is best for 9 At <^ Iwhen only one event is captured in At. The left 



16 L. Warszawski and A. Melatos 



Quantity 
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IQl^ 
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rad s-i 




10-2 




TVv 


1018 




a 


10-9 





Table 1. Fiducial pulsar parameters. 



panel of Fig. |8] graphs the PDF of n {solid black curve); a power-law fit is graphed as a dotted curve. In this example, with 
Ark ~ 1 (A is a function of 7, and hence t), p{n) covers two decades. In simulations not shown here, we investigated how p{n) 
changes with the strength of the knock-on parameter, whilst keeping Ep constant. We find that for A7/A70 < 10^, p{n) is 
sharply peaked (as in the top right panels of Fig. [3] and 2]). For 10^ < A7/A70 < 10®, p{n) extends over several decades, and 
is well described by a power law, as in the example plotted in Fig. [7] and [S] For A7/A70 > 10®, Ark exceeds unity, and hence 
every trigger event results in a system-spanning avalanche (ie. all vortices unpin). In this regime, the power law changes into 
a hump at large n. 

The fact that the system is sensitive to Ark, rather than self-organising to always produce power-law-distributed glitches, 
reminds us that the primary trigger events can decelerate the superfiuid quickly enough to keep pace with the crust; avalanches 
are not essential and only occur when A7 is large. As demonstrated in Fig. [71 however, if A7 is too large, system-spanning 
avalanches are inevitable as soon as a trigger event occurs, leading to the overpopulation of large glitch events. This behaviour 
is quantified in Sec. 18.31 Figure [9] shows that the power-law index, a, decreases for increasing A7, indicating that larger glitches 
become more probable as the knock-on mechanism strengthens. 

Waiting times as a function of glitch epoch and the waiting-time PDF are graphed in the centre panels of Fig. [T] The 
PDF spikes at small r, and there is significant statistical weight in the tail, extending out to r = 11. Excising the initial 
spike, and considering only r < 4, we fit an exponential {dotted curve) to the entire data set. p(r) = 6'(7oq) exp[— 0(7oq)r] 
is not a good fit to the data. To understand this, we consider the PDF of 7 in the bottom right panel of Fig. [7] We find 
7oq = 0.121 -C 6.95 = {7), so ^eq significantly overestimates the unpinning rate. This is not surprising since Eq. (]40p defines 
a lag that ensures that the trigger unpinning rate is enough to decelerate the superfiuid with the crust. Large jumps in the 
time series of 7 (A7/7 ~ 1, bottom left panel), corresponding to abrupt flc jumps seen in the top left panel, explain why p{'y) 
is broad (a range of 7.11 compared to 0.4 in the simple, no-knock-on case). 

The time-averaged PDF of occupied pinning energies, g{Ep), is graphed in Fig. [8] The bias towards weak pinning energies 
is in stark contrast to the no-knock-on case described in Fig.[6]in which g{Ep) is excavated. This is because the scale invariance 
promoted by 7 variability outweighs the escavation effect. 



8.3 Knock-on as a branching process 

From the Monte-Carlo simulations, we conclude that knock-on is effective in catalysing avalanche events. However, a power- 
law event size distribution only arises when the strength of the knock-on mechanism is tuned so that the boosted unpinning 
rate after an initial trigger is just sufficient to raise the probability of a subsequent unpinning to unity. This is the condition 
Ark = 1 referred to above. In this section we make this result precise using the theoretical framework of branching processes. 

From Sec.[8l we know that each trigger unpins, on average, A(i5p)rk vortices. Note that rk is the time-scale for a single 
branch in the avalanche, not the duration of the entire avalanche. The number of branches is typically the logarithm of the 
total number of vortices that unpin in the avalanche, and hence the avalanche typically lasts for < lO'^rk. 



The trigger and subsequent knock-on events form the bas is of a Gallon- Wat son-Bienayme branching process, which 
generates events that follow a generalise d Poisson di stribution ([Dobson et al.l 120071 ). The fate of familial lineages originally 
spurred the study of branching processes (|Harrislll989l l. More recently, b ranching processes have b een used to describe natural 
phenomena such a s cascades of splitting particles and p olymerisation i Tobita Sz Hamielec 1989l 'l. stati stical studies of DNA 
l Avise et al. 1984 ). the formation of random fractal sets ( Jagers 20051 ). and models of electricity grids ( Carreras et al. 2002). 

Consider a composite branching process in which the triggers follow a Poisson distribution, with mean number of triggers 
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^Ttrig in time Ttrig- Each trigger and subsequent knock-on unpinning causes, o n average, a further Ark 1) vortices to unpin. 
The event size distribution, p(n; rtrig), for the composite process is given by (jConsul fc Jainll 19731 : lOobson et allboOTf ) 



p(n; Ttrig) = ^Ttrig (nArk + druig)" 



(31) 



where n is the total number of vortices that unpin in a single avalanche. As discussed in SecO after each avalanche, the lag 
between the crust and superfluid adjusts, reducing the unpinning rate everywhere. Therefore, there is no unique waiting time 
for the trigger events, as the rate varies on the same time-scale that the triggers recur. Equation pill then implies a mean 
event size and mean unpinning rate if one waits for a time (which picks up one trigger on average) of 

M = ^A-. (32) 



1- Ark 



and 

d{n) 
dt 



(33) 



1-Ark 
respectively. 

The shape of p(n) depends strongly on A. For Ark ~ 1, which corresponds to each unpinning in an avalanche catalysing 
one other on average, p{n) exhibits a power-law tail with exponent —3/2. To see this mathematically, we use Stirling's 
/2nn (n/e)", to obtain 

„n-nATk 

1+ ■ ■ 



approximation n! 
lim p(n) 

n — ^oo 

p(n) 



A" 



nArk 



-3/2 



lim 



2tt 



(34) 
(35) 



For a state-independent system, the condition Ark = 1 leading to a power-law size distribution requires a high level of fine 
tuning. Even with the self-organising influence of torque feedback, fine tuning, in the sense of Ark ~ 1, does not arise naturally. 
Torque feedback ensures d{n) /dt is the correct value to accommodate the electromagnetic spin-down torque imposed externally, 
but it does nothing to push Ark towards unity. In other words, vortex unpinning in non-power-law-distributed avalanches is 
sufficient to decelerate the superfiuid at the same rate as the crust on average (see the top left panel of Fig. [3] for example), 
even without Ark ~ 1- 

The fine tuning condition can be re-expressed as 

PEp = -(7eq -A-y)-^ In To. (36) 

Equation (|36|l can variously be interpreted as a condition on pinning strength or stellar temperature, neither of which are 
regulated by vortex unpinning. That is, /3 is determined independently by stellar cooling, 7eq results from the mean torque, 
A7 is a direct result of quantum mechanics and the pinning landscape, and Fq is also quantum mechanical. For a typical 
pulsar, with physical parameters given in Table[l]and e = 1/2, we obtain A70 ~ 10^'* and A <C 1, which does not satisfy the 
stated condition for an avalanching system. 

An alternative collective mechanism is the accumulation of sound waves from many unpinning events, which vibrate the 
superfluid and increase the unpinning rate by raising the effective temperature, The estimate o f the e nergy emitted 

in acoustic radiation by a vortex (A_Esound = 10~^^ Jm~^) provided at the end of Warszawski et al. ( 2012h suggests that 
single unpinning events do not output sufficient acoustic energy to noticeably change the system temperature. However, the 
cumulative effect of many vortices unpinning may be enough to significantly alter the unpinning rate; of course, we are still 
obliged to identify a large-scale unpinning trigger. We defer this investigation to further work. 



9 FUTURE IMPROVEMENTS 

The framework presented in this paper can be extended in many profitable directions, two of which are outlined briefiy here. 



9.1 Towards an analytic description: master equation 



In a state-dependent system, we are interested in the probability of occupying a g iven state. The rel e vant state variable in 
the superfluid-crust system is the differential rotation, parametrised by 7. Following Daly fc Porporatol ( 2007 ). we employ the 
Chapman-Kolmogorov forward equation to track the evolution of the probability p{'y,t)d^ of being in a state 7 at time t: 



9p(7,i) _ d 



dt 



[7op(7, t)] - 9(7, t)p(7, t) + / d(A7)e(7 - A7, f)p(7 - A7, t)h{A'y, 7 - A7) 



(37) 
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The first term on the right-hand side of Eq. (|37[l describes deterministic spin down in response to the electromagnetic torque; 
we have 70 — Nc/{IcAQci)- The second term is the rate at which the system jumps away from the state 7; 0(7, t) is the 
unpinning rate at a given 7. The third term accounts for the rate at which the system jumps into the state 7 via a jump of 
size A7; here, h{A'y, 7 — A7) is the probability density of an event of size A7, when the system is in a state 7 — A7. 

For a vortex pinned with energy Ep, the unpinning rate is 6(7) ~ Foe^^^'''' . In the simplest case, discussed in Sec. 21 each 
event consists of the unpinning, outward motion, and repinning of a single vortex. There is no knock-on and hence no vortex 
avalanches. The jump size distribution in this case is h{Ay, 7 — A7) — S{A"f — A71), where A71 is the change in 7 due to the 
unpinning of a single vortex, and S{-) is the Dirac delta function. More complicated jump distribu tions can be incorporate d 
into Eq. (|37|) . via /i(A7, 7 — A7). For example, rainfall events are modelled by a gamma distribution j Dalv fc Porporatollioiol ). 
and hyperexponentials are used in queueing and communication problems ( Feldmann fc Whitt 1998h . Equation (|31|l gives the 
jump distribution appropriate to a branching process. 

We can calculate the moments of p(7). Assuming a constant A7 = A71 [i.e. h{A'y, 7 — A7) = 5{A^ — A71)], multiplying 
Eq. p7|) by 7, and integrating by parts with respect to 7, we obtain 

^ = -7o + iVroA7ir,(/3i5p,f) , (38) 
with ri(l3Ep,t) = d'ye-'^'^^''p{'y,t) = {e"''^p^). In the steady state, Eq. ^ gives 



viPEp) 



70 



iVroA7i 



(39) 



Since the right hand side of Eq. (|39[) is independent of PEp, rj must be too. Rearranging Eq. (f39 
g-^Sp(7> j-g^ good approximation if p(7,t) is sharply peaked], we arrive at 



and assuming (e 



-/9Bp7 



(7) 



1 



In 



70 
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- 7oq 



(40) 



which agrees with Eq. (|17|l . 

Solvi ng Eg. (I37II self-c onsistently is an ongoing avenue of investigation. We take inspiration from studies of solar flare 
statistics ( Wheatlandl 20091 ) , in which a Monte-Carlo algorithm, based on a power-law jump size distribution [/i( A7, 7 — A7)] , 
is employed. Analytic solutions for gamma-function and hyperexponential jump size distributions have also been described 
bv balv fc Porporatol (|201ol '). Heavier tails in h{A"f 7 — A7) lead to larger excursions from steady spin down, since there are 
more large events that push the system away from equilibrium. 



9.2 The real pulsar geometry 



This model does not include the important consideration of the radial position of the vortex, which is relevant to several aspects 
of the model. Firstly, the Magnus force experienced per unit vortex length depends directly on velocity through Eq[T] with the 
direct corollary that vortices at different radii experience a different Magnus force, even under the same lag conditions between 
the superfulid and neutron star crust. In the model presented here, we make the ansatz that glitch-relevant pinning occurs in 
only a thin spherical shell in the outer kilometer of the Neutron star, and thus have treated all vortices as having a common 
radius. Inclusion of a radial component to such a statistical model would require a corresponding probability distribution of 
radial position, whi ch would in turn demand an understanding of the radial dependence of pinning strength, as discussed by 
Haskell et all J2OI2I ). One then expects a fluctua -ting history of pile- up at different radii. Such spa. tial inhomogeneity opens up 
a new channel to drive sandpile-like avalanches ( Melatos et al. 20081 : Warszawski fc Melatos 2008) , supplem enting the pile-up 
that occurs in energy space through the coherent-noise process encoded in a(E.t) jWarszawski et al.ll2009l '). 

The second important effect of radial position is on the length of the vortex and fraction thereof that is intercepted by 
the crustal pinning region. The greater the fraction of the vortex submerged in the crust, the greater t he energy with whic h 
the vortex is pinned, resulting in an increasing function of pinning energy with radial vortex position ( Haskell et al. 20121 ). 
This separate effect would also need to be accounted for in the pinning energy distribution. 

It should be noted that the importance of tracking radial position is somewhat lessened if vortices form a turbulent tangle, 
driven by relative flow between the visco us and inviscid superfluid compon ents, as in the Donnelly-Glaberson instability 



( Peralta et ahlbooel : lAndersson et al.ll2007l ) or its pinning-mediated analogue (jLinkl l l2012al lbl). In such a scenario, small vortex 



loops in the tangle behave more like point vortices and the homogeneous statistical mo del presented in this paper is a better 
appr oximation than otherwise. However, rotation overall polarises the vortex tangle (jMongiovi et al. 2007 : Tsubota et al 



2012I ). so there is a residual geometric effect, which again would be captured by making g{E,t) depend on cylindrical radius 
in a more complete model. 
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10 SUMMARY 



In this paper we present the foundations for a statistical model of pulsar glitches. The model contains three novel features: (1) 
treatment of vortex unpinning as a Poisson process with variable rate governed by crust-superfluid lag; (2) pinning sites with 
multiple energies, which lead to excavation of the energy distribution; and (3) knock-on, which catalyses vortex avalanches. 
However, the model is incomplete because it requires fine tuning to generate scale-invariant avalanches with a power-law size 
distribution. 

In Sec. [4l we demonstrate that torque feedback between the crust and superfluid allows the angular velocity lag between 
the two components to fluctuate around a constant value. An average-rate calculation, which assumes that unpinning matches 
spin down on average, and a more comprehensive description in terms of a Chapman-Kolmogorov equation fSec. 19.1] ). predict 
the same equilibrium lag, which agrees reasonably well with results from asynchronous Monte Carlo simulations with a single 
pinning energy. Importantly, however, torque feedback alone is insufficient to catalyse a broad range of glitch sizes over many 
decades, even with a broad range of pinning energies. This result is significant, as it means that correlated vortex unpinning 
is essential to explain glitch data. 

In Sec. ISl we introduce a parametrisation of knock-on between unpinned and still-pinned vortices, which makes it easier 
for additional vortices t o unpin following an init ial unpinning. We remain agnostic about the precise origin of knock-on. Gross- 
Pitaevskii simulations ( Warszawski et al. 2012h point to (1) a vortex proximity effect, which results from an increase in the 
superfluid velocity at the site of a pinned vortex, when an unpinned vortex passes nearby; and (2) dislodgement by acoustic 
pulses emitted whenever a vortex repins, which vibrate the superfluid. Back-of-the-envelope calculations of the magnitude 
suggest that, for predicted pinning strengths and temperatures in a pulsar, neither mechanism is sufficiently strong to catalyse 
large-scale unpinning events, (i.e. to give Ark = 1). The level of fine-tuning in a model which includes knock-on is quantified. 
If kno ck-on is too strong , most glitches are system-spanning; too weak and single unpinning events dominate, akin to vortex 
creep jAlpar et al.ll 19841 ) . 

In closing, we emphasise that the Monte-Carlo simulations described in Sec. 18. II include all the new physics described in 
Sec. 12.21 whereas the master equation (Sec. 19. 1|) awaits further work. We cannot explicitly extract the fine-tuning conditions 
from the automaton. Hence a solution to Eq. (|37[) is worth pursuing. 

The lack of spatial structure in our model is an obvious deficiency. Cheng et al. ( 19881 '): Melatos et al. ( 20081 ') and many 
others have cited a richly-connected network of capacitive regions of pinned vortices as the physical origin of broadly distributed 
glitch sizes. In this way, even when Ark ~ 1, the entire system of vortices is not 'available' to the avalanche, since the knock- 
on mechanism may peter out once a local 'stress reservoir' is exhausted. We propose two ways to do this: (1) combine the 
automaton described in Sec. 18.11 with an automaton [like the one studied in IWarszawski fc MelatosI l|2008l )] that explicitly 
models nearest-neighbour interactions on a grid of pinning sites; or (2), make the assumption that strong pinning sites are 
somehow clustered (at grain boundaries in the crust, for example). We then speculate that the vortex density, and hence the 
distance of approach of an unpi nned vortex to nearby p inned vortices, is generally higher in stronger pinning regions (some 
evidence of this is described in Warszawski et al. ( 20121 )1. so that A7 is an increasing function of Ep. Broadly distributed 
glitch sizes observed in preliminary simulations that include this dependence (not presented here) are encouraging. We will 
pursue this line of thinking in future work. 
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APPENDIX A: SUPERFLUID ANGULAR MOMENTUM AND VORTEX MOTION 

Consider the angular momentum of a superfluid, projected onto the rotation axis, Lz, containing a smooth, continuous 
distribution of vortices, 

Lz = 2iTRpf drr\vxvs\, (Al) 
Jo 

where p is the superfluid density, is the superfluid velocity, and we assume that the system is uniform in the direction of 
the rotation axis, and circularly symmetric in the transverse plane. We can rewrite |vs| in terms of the number of vortices 
N{r) at radii up to r, 

_ KiV(r) 
27rr ' 

such that Eq. ([Aip becomes 

L.^2nRprdr'-^. (A3) 
Jo 27r 

Assume now that AN vortices move from radius Tii to R2, with R2 — Ri = aR, while all other vortices remain stationary. 
L2 then becomes 



(A2) 



Lz — ALz = Rtip 



( fRi rR2 rR "1 

<^ / drrN(r)+ / dr r [N (r) - AN] + / drrN{r)\ . (A4) 

iJo JRi J Rn J 



For Ri « R, we get 

ALz ^ -pnaANR^ . (A5) 



APPENDIX B: RELAXATION TO EQUILIBRIUM 



In this appendix, we solve (|2ip explicitly for g{Ep,t) given <^(-Ep) as in Eq. (|20|) . taking Fq and Ep — Eq to be constants 
This latter assumption reflects the observational fact that the glitch size and waiting-time distributions appe ar to converge 



to a quasistationary state over years, i.e. over an interval much shorter than the spin-down time-scale Q,/Q (jMelatos et al 



20081 ). The reason for solving analytically is to find explicitly the time-scales over which relaxation to equilibrium occurs. We 



Laplace transform Eq. (|21|l . defining g{Ep,s) — dt' e g{Ep,t'), and rearrange to obtain 
g(Ep, s) = — ;^ ^-f^ H 



s + Foe-'3^p s-|-Foe-'3^p 
Integrating (|B1[) over all Ep, and noting = dE' g{E' , s), we find 



Foa;(s) 



Equation (|B2I 
function (l20l). 



Foa;(s) 

and hence 
g{Ep, s) = 



dE 



, g{E',t = 0) 







s -HFoe- 



dE' ■ 



4>{E') 



(Bl) 



(B2) 



s + Foe 

is general. Let us evaluate it for the reasonable special case g[E, t = 0) = <t>{E), with cj>[E) given by the top-hat 
One then finds 



2PEo 



In 



PH{Ep)H{2Eo - Ep) 
(,s- + Foe-''-p)ln(^2!f|^) 



(B3) 



(B4) 



We recover the initial distribution (|20|l as expected from the initial value theorem g{E, t — 
recover the stationary distribution (|26p from the final value theorem, which says that g{E,t 
all the poles of sg{E,s) lie in the left half plane. We can invert the Laplace transform (|B4|) by evaluating the Bromwich 
integral 



= lims^ 
00) — lims 



o sg{E, s). We also 
+osg{E,s) because 
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1 fr+ioo 

9iEp,t) = — dse''g{E^,s) , (B5) 

^7^* J~l-ioo 

where s is now a complex number and the integral is along the line Re{s) = 7, where 7 is taken to the right of all the poles 
of e'^^g{E, s) in the complex plane. 



